function gradient = quat2rotmatGradient(quat)

qx = quat(1); qy = quat(2); qz = quat(3); qw = quat(4); 
gradient = zeros(3, 3, 4);
                 
gradient(:, :, 1) = [ 0     2*qy  2*qz; 
                      2*qy -4*qx -2*qw;
                      2*qz  2*qw -4*qx];
                 
gradient(:, :, 2) = [-4*qy  2*qx  2*qw; 
                      2*qx  0    -2*qz;
                     -2*qw  2*qz -4*qy];
                 
gradient(:, :, 3) = [-4*qz -2*qw  2*qx; 
                      2*qw -4*qz  2*qy;
                      2*qx  2*qy  0];
                  
gradient(:, :, 4) = [ 0    -2*qz  2*qy; 
                      2*qz  0    -2*qx;
                     -2*qy  2*qx  0];
